function DecodingProbability_v3()

% Finite precision sensitive :/

clear all
close all
clc

% Working example!
    % n1=12;
    % k1=10;
    % r=k1;
    % q=2;
    % dm=DifferentMatrices(n1,k1,q)
    % test_val=MatricesWithRank(n1,k1,r,q)/dm
    % % 1-test_val

%% Quick testing! Match simulation

% Gamma1 = 0.5
% Gamma2 = 0.25
% Gamma2 = 0.25
% Layer1 = 20
% Layer2 = 40
% Layer3 = 60

% Configure parameters
T=100; % nb of transmitted packets
q=2; % field cardinality
k1=10; %
k2=10; %
k3=10; %
K1=k1;
K2=k1+k2;
K3=k1+k2+k3;
g1=0.5; %
g2=0.25; %
g3=0.25; %
% n1= %
% n2= %
% n3= %

% Preparing decoding probability vectors
p_decode_l1=zeros(T,1);
p_decode_l2=zeros(T,1);
p_decode_l3=zeros(T,1);




lucky_vector=zeros(T,1);

for w=k1:T
    for q=0:9
        lucky_vector(w)=lucky_vector(w)+ksuccesntrials(w,q,g1);
        lucky_vector(w)=lucky_vector(w);
    end
end

1-lucky_vector

% Transmitting layer 1
r1=K1;
% p_decode_l1=zeros(T,1);
for g=1:T
%     if g*g1>=k1
        p_decode_l1(g)=ProbMatricesWithRank(g,k1,r1,q);
%     end
end
p_decode_l1
p_decode_l1=p_decode_l1.*(1-lucky_vector)
1-lucky_vector
% value

% Transmitting layer 2
% r2=K2;
% for g=k2:T
%     p_decode_l2(g)=ProbMatricesWithRank(g,K2-r1,r2-r1,q);
% end
% 
% p_decode_l2=p_decode_l2.*p_decode_l1;

% Transmitting layer 2
% r3=K3;
% for g=k2:T
%     p_decode_l3(g)=ProbMatricesWithRank(g,K3-r2,r3-r2,q);
% end
% 
% p_decode_l3=p_decode_l3.*p_decode_l2;



% Transmitting layer 2
% k2=20;
% r=20;
% p_decode_l2=zeros(60,1);
% for g=k2:60
%     p_decode_l2(g)=ProbMatricesWithRank(g,40-20,40-20,q)*ProbMatricesWithRank(g,k1,r,q);
% end

% p_decode_l2=p_decode_l2.*p_decode_l1
% p_decode_l2






% Quick plotting
figure(1)
hold('on')
plot(1:T,p_decode_l1)
% plot(1:T,p_decode_l2)
% plot(1:T,p_decode_l3)
hold('off')

grid('on')
pbaspect([2.5 1 1])
set(gca,'XTick',0:2:3000)
xlim([0 60])
ylim([0 1])

% Save plot
print(gcf,'uep_ew_analytic.eps')




end


% Should work
function PMWR = ProbMatricesWithRank(m,n,r,q)

% Get first set of gaussian coefficients
gc=gausscoeffs(n,r,q);

% Calculate "sum"
val=0;
for k=0:r
    % This should be the one!
    val=val+((-1)^(r-k)*gausscoeffs(r,k,q)*q^(m*k+binomcoeffs(r-k,2)-n*m));
end

% Return probability of matrix 'm'x'n' with rank 'r'
PMWR=gc*val;

end

% Deprecated (Numerical fix)
% Should work (No need to test)
% function DM = DifferentMatrices(m,n,q)
% DM=q^(m*n);
% end

% Should work (Tested! see bottom)
function GC = gausscoeffs(m,r,q)
if r==0
    % disp('r = 0 in gauss coeffs')
    GC=1;
elseif r>0
    % disp('r > 0 in gauss coeffs')
    
    % Calculate numerator
    num=1;
    for w=m:-1:m-r+1
        num=num*(q^w-1);
    end
    
    % Calculate denominator
    denom=1;
    for w=r:-1:1
        denom=denom*(q^w-1);
    end
    
    % Calculate gaussian coefficient
    GC=num/denom;
    
elseif r<0
    disp('r < 0 error in gausscoeffs!!!')
end


end

% Not tested!!! ffs!
function bc = binomcoeffs(a,k)
% As on page 123 in "A course in combinatorics"

tmp_vector=ones(2,1);
tmp_index=1;

for w=0:-1:-k+1
    tmp_vector(tmp_index)=(a+w); 
    tmp_index=tmp_index+1;
end

num=prod(tmp_vector);

denom=factorial(k);
bc=num/denom;

end






function p = ksuccesntrials(n,k,p)
p=nchoosek(n,k)*p^k*(1-p)^(n-k);
end

















%% Testing Gaussian Coefficient generater

% Correct values are from: http://mathworld.wolfram.com/q-BinomialCoefficient.html
%
% m=2;
% r=1;
% q=2;
%
% if gausscoeffs(m,r,q) == (1+q)
%     disp('gauss test 1 succes!')
% end
%
% m=3;
% r=1;
% q=2;
%
% if gausscoeffs(m,r,q) == (1+q+q^2)
%     disp('gauss test 2 succes!')
% end








%%  !!!!!!!!!!!!!!!!!junk!!!!!!!!!!!!!!!!111

%val=val+((-1)^(r-k)*gausscoeffs(r,k,q)*q^(n*k+nchoosek(r-k,2)));
% NOTE: We must not take nchoosek(0,2)!
% Are we making a mistake? The book says we should?
% Is this right?
% val=val+((-1)^(r-k)*gausscoeffs(r,k,q)*q^(n*k+nchoosek(r-k,2)));
